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ABSTRACT 

The  Rapid  Solution  of  the  Laplace  Equation  on  Regions 

with  Fractal  Boundaries 


Jin  Hong  Ma 
Yale  University 
1992 


Interest  in  the  numerical  solution  of  the  Laplace  equation  on  regions  with  fractal 
boundaries  arises  both  in  mathematics  and  physics.  In  mathematics,  examples  in¬ 
clude  harmonic  measure  of  fractals,  complex  iteration  theory,  and  potential  theory.  In 
physics,  examples  include  Brownian  motion,  crystallization,  electrodeposition,  viscous 
fingering,  and  diffusion-limited  aggregation.  In  a  typical  application,  the  numerical 
simulation  has  to  be  on  a  very  large  scale  involving  at  least  tens  of  thousands  of  equa¬ 
tions  with  as  many  unknowns,  in  order  to  obtain  any  meaningful  results.  Attempts 
to  use  conventional  techniques  have  encountered  insurmountable  difficulties,  due  to 
excessive  CPU  time  requirements  of  the  computations  involved.  Indeed,  conventional 
direct  algorithms  for  the  solution  of  linear  systems  require  order  0{N^)  operations  for 
the  solution  of  an  N  x  N—  problem,  while  classical  iterative  methods  require  order 
0{N'^)  operations,  with  the  constant  strongly  dependent  on  the  problem  in  question. 
In  either  case,  the  computational  expense  is  prohibitive  for  large-scale  problems.  We 
present  a  direct  algorithm  for  the  solution  of  the  Laplace  equation  on  regions  with  frac¬ 
tal  boundaries.  The  algorithm  requires  0{N)  operations  with  a  constant  dependent 
only  on  the  geometry  of  the  fractal  boundaries.  The  performance  of  the  algorithm 
is  demonstrated  by  numerical  examples,  and  applications  and  generalizations  of  the 
scheme  are  discussed. 
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Chapter  1 


Introduction 


During  the  last  decade,  the  numerical  solution  of  the  Laplace  equation  on  regions  with 
fractal  boundaries  has  been  becoming  increasingly  popular  both  in  mathematics  and 
physics.  In  mathematics,  examples  include  harmonic  measure  of  fractals,  complex 
iteration  theory,  and  potential  theory.  In  physics,  examples  include  growth  phenom¬ 
ena  such  as  crystallization,  electrodeposition,  viscous  fingering,  and  diffusion-limited 
aggregation,  where  the  harmonic  measure  governs  the  growth  of  the  fractal  surfaces 
[56].  Thus,  much  recent  work  has  been  focused  on  the  study  of  the  metric  properties 
of  harmonic  measure  on  fractals  [5],  [12],  [33],  [38],  and  [34]. 

1.1  Background 

Carleson  proved  recently  in  [12]  that  the  dimension  of  the  support  of  harmonic  mea¬ 
sure  for  any  two-dimensional  Cantor  set  is  strictly  less  than  one  .  However,  the  actual 
values  for  particular  sets  have  not  been  determined,  and  it  is  unclear  how  they  can 
be,  without  some  form  of  computer  experimentation.  In  the  behavior  of  harmonic 
measure  for  Cantor  sets  is  completely  unknown.  Thus,  several  attempts  have  been 
made  during  the  la.st  several  years  to  solve  such  problems  numerically  (see  [5]  and 
[34]). 


1 


CHAPTER  1.  INTRODUCTION 


2 


There  are  two  approaches  to  the  study  of  the  metric  properties  of  the  harmonic 
measure  on  fractals. 

1.  Viewing  the  harmonic  measure  as  the  relative  hitting  probability  at  points  on 
the  surface,  and  using  the  Monto  Carlo  method  to  conduct  computer  simulations 
on  parallel  machines  such  as  the  Connection  machine  (see  [5]). 

2.  Formulating  the  problem  as  an  integral  equation  of  the  first  kind,  and  using 
brute  force  to  solve  it  numerically. 

While  the  first  approach  has  produced  some  significant  results  (see  [5]),  the  com¬ 
putation  becomes  prohibitively  expensive  when  high  accuracy  is  desired,  due  to  the 
slow  convergence  of  the  Monto  Carlo  method  (as  is  well  known,  the  error  of  a  Monto 
Carlo  simulation  decays  like  1/\/N,  where  N  is  the  number  of  trials). 

On  the  other  hand,  the  second  approach  has  also  encountered  insurmountable 
difficulties,  due  to  excessive  CPU  time  requirements  of  the  computations  involved. 
Indeed,  in  order  to  obtain  mathematically  meaningful  results,  systems  of  linear  equa¬ 
tions  have  to  be  solved,  involving  at  least  tens  of  thousands  of  equations  with  as  many 
unknowns.  Conventional  direct  algorithms  for  the  solution  of  linear  systems  require 
order  0{N^)  operations  for  the  solution  of  an  N  x  N—  problem,  while  classical  itera¬ 
tive  methods  require  order  0{N^)  operations,  with  the  constant  strongly  dependent 
on  the  problem  in  question.  In  either  case,  the  computational  expense  is  prohibitive 
for  large-scale  problems. 

We  present  a  direct  algorithm  for  the  rapid  solution  of  the  Laplace  equation  on 
regions  with  a  certain  type  of  fractal  boundaries.  The  algorithm  requires  0{N)  op¬ 
erations  with  a  constant  dependent  only  on  the  geometrical  property  of  the  fractal 
boundaries,  where  iV  is  the  number  of  elements  in  the  discretization  of  the  fractal. 
And  the  evaluation  of  the  potential  at  any  point  requires  0(log(.V))  operations. 
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1.2  Genealogy  of  Fast  Algorithms 

The  algorithm  of  this  thesis  is  closely  related  to  other  analysis-based  fast  algorithms. 
Among  them,  perhaps  the  Fast  Multipole  Method  (FMM)  [25]  is  the  best  known.  It 
provides  a  fast  scheme  to  evaluate  gravitational  and  electrostatic  potentials  involving 
a  large  number  of  particles.  In  [8],  wavelets  are  used  to  generalize  the  FMM  to  a 
variety  cf  integral  operators.  In  [4],  -wavelets  are  combined  with  the  Schultz  algorithm 
to  solve  integral  equations  of  the  second  kind.  Several  other  “faist”  schemes  have  been 
constructed,  such  as  the  algorithm  of  [27]  for  the  rapid  evaluation  of  Gauss  transforms. 
In  [48],  a  fast  direct  solver  is  developed  for  the  solution  of  integral  equations  in  one 
dimension.  However,  the  latter  algorithm  substantially  exploits  the  one-dimensional 
geometry.  While  the  algorithms  of  [25],  [8],  [4],  [3],  and  [27]  are  essentially  fast 
algorithms  for  applying  a  matrix  to  a  vector,  our  algorithm  can  be  viewed  as  a  fast 
direct  inversion  scheme  for  the  some  of  matrices  of  the  type  produced  by  the  FMM. 

1.3  Outline  of  the  Dissertation 

The  direct  algorithm  presented  in  this  thesis  exploits  the  fact  that  far-held  interactions 
are  of  low  rank  for  any  given  precision,  and  low  rank  operators  can  be  recursively 
compressed  without  actually  generating  them. 

We  begin  with  the  dehnition  of  the  problems  to  be  addressed  in  Chapter  2.  In 
Chapter  3,  -we  summarize  certain  mathematical  and  numerical  facts  to  be  used  in  this 
thesis.  In  Chapter  4,  we  establish  the  principal  analytical  tool  of  this  thesis  that  ranks 
of  far-held  interactions  are  hnite  to  any  given  precision.  In  Chapter  5,  we  develop 
the  mathematical  apparatus  used  to  construct  the  fast  algorithm  by  borrowing  termi¬ 
nology  from  the  standard  scattering  theory  for  the  Helmholtz  equation.  In  Chapter 
6,  we  present  the  description  of  the  feist  algorithm,  and  in  Chapter  7,  we  illustrate 
the  performance  of  the  algorithm  by  numerical  examples.  Finally,  in  Chapter  S.  we 
outline  some  applications  and  generalizations. 


Chapter  2 

Statement  of  the  Problems 


As  is  well-known,  the  governing  equation  for  potential  problems  is  the  Laplace  equa¬ 
tion 


d^u  d'^u 


(2.1) 


Functions  which  satisfy  (2.1)  are  referred  to  as  harmonic  functions. 

In  this  chapter,  we  define  the  problems  to  be  addressed,  namely,  the  boundary 
value  problems  for  the  Laplace  equation  on  regions  with  fractals  of  Cantor  type  as 
the  boundaries. 


2.1  Fractal  Boundaries 

A  fractal  of  (  ator  type  is  a  classical  example  of  fractals  (see,  for  example,  [7],  [16], 
and  [39]),  which  can  be  generated  recursively  by  dividing  a  given  region  into  four 
corner  regions  (boxes)  with  a  ratio  of  sides  as  a  parameter. 

Given  a  real  number  a  (  0  <  a  <  |  ),  we  define  a  sequence  of  sets  as  follows  (See 
Figures  2.1,  2.z,  and  2.3): 

Cq  =  {  the  unit  square  },  ('“•2) 

=  {4  corner  boxes  with  as  their  sizes  }.  (2.3) 
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Figure  2.1:  Set  Cq 


Figure  2.2:  Set  Cf 


□D 

□  □  □□ 


Figure  2.3:  Set  C2 
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=  {4*  corner  boxes  with  a‘  as  their  sizes  },  (2-4) 

where  /  is  an  integer. 

The  Cantor  set  associated  with  the  ratio  a  is  the  limit  of  the  sequence  of  sets 
{C“},  which  decreases  monotonically: 

The  limit  is  defined  as  the  intersection  of  the  sets  {C“} 

C"  =  fl  Cf .  (2.5) 

/e{0.1,2,..,oo} 

For  a  given  /  >  1,  we  will  define  the  /-th  approximation  to  the  Cantor  set  (asso¬ 
ciated  with  the  ratio  a)  as  a  set 

At  =  {  the  centers  of  all  boxes  in  Cf  }.  (2.6) 

We  will  refer  to  the  4^  boxes  Cf  generated  during  the  /  —  th  step  of  the  above 
process  as  level  I  boxes.  Thus,  there  is  one  box  on  level  0,  and  it  coincides  with  the 
unite  square.  The  level  /  -|-  1  is  obtained  from  the  level  I  by  subdividing  each  box  on 
the  level  /  into  4  corner  boxes  (see  Figure  2.2). 

We  will  also  impose  a  tree  structure  on  the  hierarchical  structure  of  C®,  so  that  if 
ibox  is  a  fixed  box  at  level  1,  the  four  boxes  at  level  /  -)-  1  obtained  by  subdivision  of 
ibox  are  considered  its  children,  while  the  four  child  boxes  are  considered  neighbors. 


2.2  The  Laplace  Equation 

Let  C“  denote  the  Cantor  set  associated  with  the  given  ratio  a  (see  Figure  2.4).  NIe 
will  consider  the  following  exterior  Dirichlet  problem  for  the  Laplace  equation 

Au  =0  for  X  E  \  C“, 
u\c‘^  =  /. 
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Figure  2.4:  Cantor  Set  C“  -  A  Fractal 
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To  insure  the  uniqueness  of  the  solution  of  problem  (2.7),  the  far-held  condition 

lim  u(x)  —  rlnr^  ^  ^  =0  (2-8) 

r—oo 

is  normally  imposed.  In  the  above  formula,  =  (a:  —  xq)^  A  (y  —  J/o)^  with  point 
X  =  (x,y)  €  and  an  arbitrary  fixed  point  xq  =  (io,2/o)  G  C’“- 

Remark  2.1  .4s  is  well-known,  for  any  harmonic  function  u  :  R?  ^  R} ,  there  exists 
two  functions,  and  -tf, 

OO 

k——oo 

00 

-f- 7lnr, 

k=—<x> 

such  that  u  =  (f  tl). 

The  far-field  condition  (2.8)  excludes  the  constant  term  in  the  expansion  of  ip 
while  allowing  the  logarithm  term. 

The  proof  of  the  following  theorem  can  be  found,  for  example,  in  [55]. 

Theorem  2.1  The  boundary  value  problem  (2.7)  with  the  far-field  condition  (2.8)  is 
a  well-posed  problem. 

As  is  well-known  (see  [15],  for  example),  the  boundary  value  problem  (2.7)  with 
the  fax-field  condition  (2.8)  can  be  formulated  as  an  integral  equation  of  the  first  kind 
by  representing  the  solution  as  the  logarithmic  potential  of  the  charge  distribution  on 
the  boundary  C“.  The  charge  distribution  is  a  Borel  measure  on  the  Cantor  set 
(see  [15]).  Denoting  by  a  the  charge  distribution  over  the  boundary  C“,  we  obtain 
the  integral  equation 

f  In  |x  —  f  |d<T(f)  = /(i)  (2.9) 

JC’^ 

with  I  €  C“,  where  the  integration  is  in  the  sense  of  Borel  measure. 

Remark  2.2  As  is  well-known,  the  Lebesgue  measure  of  a  Cantor  set  is  zero. 
However,  the  equation  (2.9)  is  mathematically  sound  due  to  the  fact  that  any  Cantor 
set  C“  possesses  a  positive  capacity  (for  detail,  see  [1.5]  or  [55]) 
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For  a  given  integer  Z,  >  1,  suppose  that  N  =  4^,  and  =  {2,  |  z  =  1, 2, . . . ,  iV} 
is  the  L-th  approximation  to  the  Cantor  set  C“.  Then,  by  the  definition  of  Borel 
integration,  the  integral  equation  (2.9)  is  discretized  as  a  linear  system  of  equations 

A(7  =  b  with  AeR^^^  and  b  e  R^ ,  (2.10) 

where 

Aij  =  log  I  2,'  —  2j  1  for  i  ^  j,  and  Aij  =  /  •  log  a  +  6  for  i  =  j.  (2.11) 

In  formula  (2.11),  6  is  a  constant  dependent  on  the  ratio  a  and  the  total  number  of 
levels  L  in  the  approximation  of  Cantor  set  C“.  The  actual  choice  of  the  constant  8 
will  be  discussed  in  Chapter  7. 

Given  a  measure  fi  over  C“,  we  define  an  integral  operator  P^^  :  R^  R}  by  the 
formula 

P^(x)=/  ln|x-t|d/z(t).  (2.12) 

JC‘^ 

Lemma  2.1  can  very  easily  verified;  and  it  states  that  with  any  non-negative 
measure  fx  over  (7“,  the  function  P^{x)  is  negative  for  x  €  C“.  Theorem  2.2  is  the 
immediate  consequence  of  Lemma  2.1. 

Lemma  2.1  If  n  is  a  non-negative  measure  overC°^,  then 

P^{x)<0  forxeC\  (2.13) 

Theorem  2.2  There  exists  an  integer  Nq  such  that  for  any  integral  N  >  No,  the 
coefficient  matrix  A  in  (2.10)  is  negative  definite. 

Remark  2.3  Experiments  show  that  the  integer  No  in  Theorem  2.2  can  be  as  small 
as  No  =  64  (see  Chapter  7). 


Chapter  3 


Mathematical  and  Numerical 
Preliminaries 


In  this  chapter,  we  summarize  certain  well-known  mathematical  and  numerical  facts 
to  be  used  in  the  rest  of  this  thesis.  They  can  be  found,  for  example,  in  [10],  [13], 
and  [55]. 


3.1  Potential  Theory  for  the  Laplace  Equation 

3.1.1  Green’s  Function 

Definition  3.1  ^Green’s  function  in  R^)  Suppose  that  U  C  is  an  open  con¬ 
nected  set,  then  for  P  €  fl,  function  G{z,P)  ;  fl  — »  will  be  referred  to  as  the 
Green’s  function  ofQ  with  pole  (or  singularity)  at  P  if  it  satisfies  the  following  three 
conditions. 

1.  G{z,P)  is  harmonic  in  fl,  except  at  the  point  P. 

2.  If  P  ^  oo,  then 

g{z)=^G{z,P)-\og\z-P\  (3.1) 

is  harmonic  in  a  neighborhood  of  P. 


10 
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IfP  =  oc,  then 

g{z)  =  G{z,oo) -\og\z\  (3.2) 

is  harmonic  in  a  neighborhood  of  oo. 

3.  j4s  2  tends  to  any  point  on  the  boundary  ofD,  then 

G(2,P)^0.  (3.3) 

Theorem  3.1  (Existence  of  Green’s  Function  in  R^)  Suppose  that  Q  C  R^  be  a 
domain  bounded  by  a  Jordan  curve  F.  Then  the  Green’s  function  G{z,P)  for  LI  exists 
for  any  P  E  Ll. 

Lemma  3.1  (Green’s  Formulaj  Let  7  :  [0,  L]  Rf  be  a  closed  Jordan  curve, 
with  image  of  7  denoted  by  F.  Suppose  that  LI  C.  Rf  is  the  interior  of  F,  so  that 
F  =  dLl.  Suppose  further  that  Nt  :  [0,F]  — >  R^  is  the  interior  normal  to  F,  function 
G:  Ll  X  Ll  R}  is  the  Green’s  function  for  LI,  and  9  €  L‘^{T).  Then  the  function 
u  :  Ll  R  defined  by  the  formula 

"  s  r 

is  harmonic,  and  ujr  =  <p- 

3.1.2  Boundary  Value  Problems  for  the  Laplace  Equation 

Suppose  that  F  C  is  a  Jordan  curve,  parameterized  by  its  length  7  :  [0,  Z]  — >  R?, 
and  Ll  is  the  region  bounded  by  F,  so  that  dLl  =  F.  Suppose  further  that  N  :  [0,  L]  — > 
is  the  interior  normal  to  F.  For  an  integrable  function  /  :  [0,  Z]  — >  we  will  be 

solving  one  of  the  following  problems. 

(A)  Interior  Dirichlet  problem 


A$(i)  =  0 
4>(x)  =  f{l~'^{x)) 


for  X  ^  Ll 
for  X  €  F 


(3.5) 
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(B)  Exterior  Neumann  problem 

A^(a;)  =  0  for 

^^(x)  =  /(7-‘(x))  for 

with  ^  satisfying  the  far-held  condition  (2.8). 

As  is  well  known,  each  of  the  above  two  problems  has  a  unique  solution  for  any 
continuous  right  hand  side  /  and  piecewise  smooth  boundary  F  (see,  for  example, 

[13]). 

3.1.3  Single  and  Double  Layer  Potentials 

Suppose  that  a  point  charge  of  unit  intensity  is  located  at  the  point  xq  G  R^.  Then, 
for  any  x  G:  with  x  ^  Xo,  the  potential  due  to  this  charge  is  described  by  the 

expression 

= -Indk  -  2;o||).  (3.7) 

The  potential  of  a  dipole  of  unit  intensity  located  at  Xq  and  oriented  at  the  direc¬ 
tion  h  G  R^  (l|/i||  =  1)  is  described  by  the  formula 

(3-8) 

||x  -  Xoll 

For  an  integrable  function  /x  :  [0,i/]  — »  the  potential  of  a  single  layer  with 

density  /x  is  given  by  the  formula 

'I'(x)  =  f  <i)^(t){x)n{t)dt,  (3.9) 

Jo 

and  the  potential  of  a  double  layer  with  the  dipole  density  /x  is  given  by  the  formula 

$(x)  =  f  (3.10) 

Jo 

3.1.4  Integral  Equations  of  the  Classical  Potential  Theory 

In  the  classical  potential  theory,  the  interior  Dirichlet  problem  (3.5)  is  solved  by 
representing  $  as  the  potential  of  a  double  layer,  and  the  exterior  Neumann  problem 


X  e  R^\n 
xgT 


(3.6) 


CHAPTER  3.  MATHEMATICAL  AND  NUMERICAL  PRELIMINARIES 


13 


(3.6)  is  solved  by  representing  ^  as  the  potential  of  a  single  layer.  The  analysis  of 
the  single  layer  and  double  layer  potentials  in  the  vicinity  of  the  boundaries  results 
in  two  integral  equations  of  the  second  kind. 

(Al)  Interior  Dirichlet  Problem 


Trn{x)  +  / 
Jo 


L  d 
dNit) 


log  il7(x)  -  'y(t)\\fi{t)dt  =  fix). 


(3.11) 


(Bl)  Exterior  Neumann  Problem 


dW{^)  lo  ^°Sll7(^)-7(0!l/^(0<^^  = /(3^)-  (3.12) 


3.2  Galerkin  Method  for  the  Solution  of  Integral 
Equations 

As  is  well-known,  the  classical  Galerkin  method  can  be  used  for  the  numerical  solution 
of  the  integral  equations  of  the  form 


fi{x)  A  j  Kix,t)fi{t)dt  =  f{x). 


(3.13) 


The  following  is  the  brief  description  of  the  Galerkin  method. 

Suppose  that  {Pi(i),  P2(^))  •  ’ '  ? -^71(3;),  •  •  •}  is  the  orthonormal  basis  in  L‘^[a,b]. 
Then  the  function  :  [0,  L]  —*  R}  defined  by  the  formula 


and  satisfying  the  condition 


/^n(x)  =  Y,QjPj{x), 
j=l 


=  0, 


(3.14) 


(3.15) 


will  be  used  to  approximate  the  solution  of  the  integral  equations  (3.13).  The  error 
function  r(x)  in  formula  (3.15)  is  defined  via  the  e.xpression 


'{x)  =  Unix)  A  K{x,t)fln{t)dt  -  fix). 


(3.16) 
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The  above  procedure  results  in  a  linear  system  Bx  =  b  defined  by  formulae 

Bij  =  f  f  K{x,t)Pi{t)Pj{t)dxdt  i  ^  ji  (3-17) 

Jo  Jo 

for  1  <  i,  j  <  n,  and  Bn  =  1, 

bi  =  f  f{x)Pi{x)dx,  (3.18) 

Jo 

for  1  <  t  <  n. 

A  detailed  discussion  of  the  Galerkin  method  for  the  solution  of  the  equations  of 
the  classical  potential  theory  can  be  found,  for  example,  in  [6],  [10],  [14],  [20],  [24], 
[32],  [52],  [53],  and  [61]. 


Chapter  4 

Fields  of  Charges 


In  this  chapter,  we  investigate  in  some  detail  the  structure  of  potential  fields  in 
and  prove  Theorem  4.2,  which  is  the  principal  analytical  tool  of  this  thesis. 

It  is  well  known  that  the  potential  (^Xo  <iue  to  a  point  charge  at  xo  €  (defined 
by  formula  (3.7))  is  harmonic  in  any  region  excluding  the  source  point  xq.  Moreover, 
for  any  harmonic  function  u  :  R^  R},  there  exists  an  analytic  function  iv  :C  —*C 
such  that  u(x,y)  =  Re{w{x,y)).  In  the  rest  of  this  thesis,  we  will  make  no  distinction 
between  points  in  R^  and  points  inC.  In  complex  terms,  the  potentials  d>Xo  and  d>Xo,ft 
defined  by  the  expressions  (3.7)  and  (3.8)  respectively,  eissume  the  form 


<t>zo{^)  =  Re{-\n{z  -  Zo)), 

and 

K.h(z)  =  fie(— ^), 

Z  -  Zo 

where  z  =  x  +  iy  and  zq  =  xo  +  iyo-  Following  the  standard  practice,  we  will  refer  to 
the  analytic  function  ln(z  —  zq)  as  the  potential  at  the  point  z  gC  due  to  a  charge 
located  at  the  point  zq. 


1-5 
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4.1  Ranks  of  Interaction  Matrices 

The  following  lemma  can  be  easily  proved  by  expanding  ln(l  —  u-’)  into  Taylor  series 
with  respect  to  u>. 

Lemma  4.1  Let  a  unit  point  charge  be  located  at  Zq-  Then  for  any  z  such  that 


l-l  >  kol, 

<?zo{2)  -  ln(2  -  2o)  =  aoln(2)  +  Tr 

(4.1) 

k=l  " 

where 

Go  =  1  awd 

(4.2) 

Furthermore, 

for  any  p  >  1, 

p 

(pzoiz)  -Ooln(2r)  -  £  3 
*=1 

(4.3) 

where 

(4.4) 

2o 


Remark  4.1  The  above  lemma  can  be  reformulated  in  a  slightly  different  way,  which 
will  be  used  later. 

Truncating  the  expansion  (4-1)  after  p  terms  (p  >  1),  we  will  denote  the  error  of 
the  truncated  expansion  by  sl^{z),  so  that 

k-1  " 


Then  for  any  p  >  I, 


Ozo{-)  = 

(4.5) 

with  the  vectors  Up  and  Vp  defined  by  the  formulae 

1  1  1 

(4.6) 

Up  =  (Inr,-,— ,•■•.— ) 

Z 

-,2 

tp  U,  ^  ,  2  •  '  P  ’■ 

(4.7) 
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Furthermore, 

(4.8) 

with  c  defined  by  (4-4)- 

Clearly,  the  truncation  error  in  (4.8)  decays  exponentially  as  the  function 

of  p.  Thus  few  terms  in  (4.1)  are  needed  to  achieve  any  given  accuracy, 

The  following  theorem  is  the  principal  analytical  tool  of  this  thesis.  It  states  that 
under  certain  conditions,  the  rank  of  the  matrix  describing  the  potential  interaction 
between  two  sets  in  is  finite  to  any  given  precision. 

We  define  the  interaction  matrix  between  the  two  sets  of  points  {xi,  ■  ■  ■  ,Xm}  and 
[Vi-,  •  •  •  iVn}  as  a  m  X  n  matrix 

••• 

<f>y2{^2)  •••  d>yn{=^2) 

0yi(^m)  '''  Pyni^'m) 

In(xi-pi)  ln(xi-?/2)  •••  ln(xi-t/„) 

ln(x2-j/i)  ln(x2-T/2)  •••  ln(x2-y„) 

ln(x„-yi)  ln(x,„-t/2)  •••  ln(x„,-i/„) 

The  following  theorem  follows  immediately  from  formulae  (4.5  )  and  (4.8). 

Theorem  4.1  Let  n  unit  point  charges  be  located  within  the  circle  \y  \  <  R  at  points 
{j/i5  J/2i  ■  ■  ■  •  J/u} •  A  >  0  6e  some  real  number,  and  {xj.  X2,  •  •  • ,  x,n)  be  another  set  of 
points  such  that  |x,|  >  (1  +  X)R  for  all  \  <  i  <  m  (see  Figure  4-1).  Then  the 
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Inequality  (4.12)  means  that  every  element  of  the  matrix  of  truncation  error  E^ 
decays  exponentially  as  the  function  of  p.  Thus  for  any  given  accuracy,  the  interaction 
matrix  of  the  two  sets  {i,}  and  can  be  decomposed  into  the  product  of  two 

matrices  of  low  rank  (  <(P+1)  )■ 
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Figure  4.2:  Well-separated  sets  in  the  plane 

For  two  sets  of  points.  {j1.x2.---.Xm}  cf  and  {j/i.  3/2.  •  •  •  ^ !/«}  we  say  that 
the  two  sets  are  well-separated  (see  Figure  4.2)  if  there  e.xist  points  xo-t/o  6(F.  and 
/?  >  0  such  that 

lx,  —  Xol  <  R  for  1  <  t  <  m, 

\yj  -  yo\  <  R  for  1  <  ;  <  n. 

i-To  -  3/ol  >  3/?. 

Theorem  4.1  implies  that  with  any  prescribed  precision,  the  interaction  matrix 
of  two  well-separated  sets  can  be  decomposed  inio  a  nrodc'-t  of  two  matrices  of  low 
rank,  the  rank  depending  only  on  the  separation  of  the  two  sets  (A  >  1).  (see  formula 
(4.12)). 


4.2  Interactions  in  Cantor  Sets 

In  this  section,  we  consider  electrostatic  interactions  within  a  Cantor  set  C“.  For  any 
given  ratio  a.  the  interactions  in  the  Cantor  set  are  of  low  rank.  The  interaction 
ranks  depend  only  on  the  ratio  a  for  generating  the  Cantor  set. 

The  following  lemma  is  obvious,  and  will  be  used  in  the  proof  of  Theorem  4.2. 

Lemma  4.2  Suppose  that  and  D2  are  two  subsets  of  Cantor  set  C“  with  ratio  a, 
set  A  is  a  child  box  of  D\.  and  B  is  a  child  box  of  D2  (see  Figure  4-3) ■  Then  the  rank 
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Figure  4.3:  Two  Subsets  of  C“  and  Their  Child  Boxes 

of  the  interaction  matrix  between  subsets  D\  and  D2  at  most  four  times  as  large  as 
that  between  boxes  A  and  B. 

Theorem  4.2  For  a  given  real  number  a(0<a<^),  and  an  integer  I  >  I, 
suppose  that  C“  ’S  the  Cantor  set  associated  with  the  ratio  a,  and  Cf  if  '  =  •  set  of 
boxes  generated  at  the  l—th  level 

=  (4-13) 

Then  for  any  given  precision,  the  rank  of  the  interaction  matrix  between  any  two  boxes 
Di  and  Dj  depends  only  the  ratio  a  for  generating  the  Cantor  set  and  does  not 
depend  on  the  sizes  of  boxes  and  the  numbers  of  points  inside  the  boxes. 

In  other  words,  the  matrix  of  interactions  between  any  two  boxes  at  any  level  of 
Cantor  set  C“  is  of  fixed  rank,  to  any  prescribed  precision. 

Proof:  Because  of  the  self-similarity  of  the  Cantor  set  C“,  it  is  sufficient  to  estimate 
the  rank  of  interactions  between  any  two  boxes  of  C“  with  /  =  1, 

C^  =  {Di,D2,Ds,D4}. 

Suppose  that  Di  and  £>2  are  not  well-separated.  Then  we  divide  each  of  them  into 
four  squares  of  the  same  size  (see  Figure  4.4).  Let  Abe  a.  square  from  the  subdivision 
of  Di,  and  S  be  a  square  from  the  subdivision  of  D2. 
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Figure  4.4:  Subdivisions  of  Two  Sets 


If  squares  A  and  B  are  well-separated,  then  the  rank  of  interaction  between  Di 
and  Dj  is  no  greater  than  4p,  where  p  is  the  interaction  rank  between  A  and  B. 

If  A  and  B  are  not  well-separated,  then  we  keep  dividing  each  region  into  four 
squares  till  the  pieces  are  well-separated. 

Due  to  Lemma  4.2  and  Theorem  4.1,  the  rank  of  the  interaction  matrix  between 
any  two  boxes  at  any  level  of  the  Cantor  set  is  bounded  by 


p  ■  4 


k 


with 


k  = 


ln(l  -2a)/(\/2-  1) 
In  a 


(4.14) 


where  p  is  the  rank  of  the  interaction  matrix  between  two  well-separated  boxes.  ■ 

Remark  4.2  Clearly,  the  estimate  (4-14)  is  an  extremely  pessimistic  one.  In  the 
following  chapter,  we  obtain  much  sharper  numerical  estimates  (see  section  5.3). 


Chapter  5 


Scattering  Theory  for  the  Laplace 
Equation 


In  this  chapter,  we  borrow  terminology  from  the  standard  scattering  theory  of  wave 
equations  for  the  design  of  fast  algorithms,  and  refer  to  the  result  as  scattering  theory 
for  the  Laplace  equation. 

To  develop  the  scattering  theory  for  the  Laplace  equation,  we  first  introduce  the 
concept  of  scattering  matrix,  and  then  present  the  merging  scheme  for  generating 
scattering  matrices  recursively. 

Throughout  this  section,  F  will  denote  a  Jordan  curve,  parameterized  by  its  length 
7  :  [0,  L\  —*  B?.  The  region  bounded  by  F  will  be  denoted  by  Q,  and  D  will  denote  a 
compact  subset  of  fl.  In  addition,  G:  B}  will  denote  the  Green’s  function  for 

domain  fl,  and  N  :  [0,  L]  —>  B?  will  denote  the  interior  normal  to  F.  For  a  compact 
set  E  C  B?,  M{E)  will  denote  the  set  of  all  non-negative  Borel  measures  on  E. 


5.1  Scattering  Matrices 

Any  function  ^  :  Vt  B  harmonic  inside  fl  and  continuous  on  will  be  referred  to 
as  incoming  potential.  As  is  well-known,  for  any  continuous  function  ip  :  T  B, 
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Figure  5.1:  Compact  set  D  in  domain  with  boundary  F 

there  exists  a  unique  function  ^  :  Q,  R  harmonic  on  fl,  and  continuous  on  Q  such 
that  $|r  =  ^p.  Therefore,  we  will  abuse  the  notation  by  referring  to  the  function 
V? :  r  -4  /2  as  the  incoming  potential. 

Suppose  that  q  €  and  cr  is  a  Borel  measure  over  D.  Given  a  function 

K  €  L^{F(?  X  R?),  the  function 

’I'(x)  =  f  K{x,t)q{t)da{t)  ior  x  £  R^\D  (5.1) 

Jd 

will  be  referred  to  as  outgoing  potential.  Similarly,  we  will  call  its  restriction 
V’  =  ’J'lr  onto  r  an  outgoing  potential.  Outside  the  domain  fl,  function  will  also 
be  referred  to  as  scattering  potential. 

Remark  5.1  Particularly,  we  are  .nterested  in  the  case  when  K{x,t)  =  In  ||x  —  t|l, 
and  q(t)  is  the  characteristic  function  of  D.  Then  the  outgoing  potential 

’P(x)  =  f  In  ||x  —  <||dcr(t)  (5.2) 

J  D 


is  a  function  harmonic  in  Rf\D,  and  satisfying  the  far  field  condition  (2.8). 
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We  define  three  operators 

L  :  L%V)  ^  L^{D), 

P  ;  MiD)  L^iR'^XD), 
S  :  M{D)  -> 

via  the  formulae 


L{p){x)  = 

(5.3) 

^ _ Q 

II 

0., 

K(x,t)q{t)d<7{t)  ioT  X  E  R^\D, 

(5.4) 

S{(t){x)  = 

1  K{x,t)q{t)da{t)  for  x  6  F. 

JD 

(5.5) 

We  will  be  considering  equations  of  the  form 

P{<^)  =  /, 

(5.6) 

with  /  €  L^{D).  A  special  case  of  equation  (5.6)  is  the  integral  equation  (2.9)  defined 
in  section  2.2,  with  D  a  Cantor  set,  K{x^t)  =  In  ||x  —  <1|,  and  q{t)  the  characteristic 
function  of  D. 

Definition  5.1  The  operator  a  :  I<^(r)  — ^  L^{T)  defined  hy  the  expression 

a  =  SP-^L  (5.7) 

will  be  referred  to  as  scattering  matrix. 

Remark  5.2  Given  an  incoming  potential  on  the  boundary  F,  the  operation  of  a 
on  ip  can  be  viewed  as  consisting  of  three  steps: 

1.  The  operator  L  constructs  a  function  f  =  Lp  .  H  —*  R  harmonic  over  the 
compact  set  D,  and  such  that  iLp)\r  =  p. 
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Figure  5.2:  Scattering  matrix  a:  (y?  — >  0 

2.  The  operator  P~^  constructs  the  solution  a  =  P~^f  of  equation  (5.6)  from  the 
harmonic  function  f  =  L^p. 

The  charge  distribution  a  will  he  referred  to  as  induced  charge  distribution. 

3.  The  operator  S  defined  by  (5.5)  constructs  an  outgoing  potential  =  Scr  on  F 
from  the  induced  charge  distribution  a . 

The  outgoing  potential  ip  will  be  referred  to  as  induced  outgoing  potential. 

Thus,  the  scattering  matrix  a  maps  an  incoming  potential  p  to  the  induced  out¬ 
going  potential  xp 

ip  =  ap.  (5.8) 

The  following  theorem,  while  interesting  in  itself,  is  not  closely  related  to  the 
purpose  of  this  thesis  since  its  proof  is  quite  involved.  We  refer  the  reader  to  [15]  and 
[55],  where  it  can  be  found  in  a  somewhat  different  form. 

Theorem  5,1  (Compactness  of  Scattering  Matrix  q)  Suppose  that  K{x,t)  = 
In]]x  —  t|],  and  q{t)  is  the  characteristic  function  of  D.  //  Z)  Pi  F  =  0,  then  the 
scattering  matrix  a  in  (5.7)  is  a  compact  operator. 
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Figure  5.3:  Disjoint  compact  sets  {A}  domain  Cl 


In  other  words,  if  D  is  strictly  inside  domain  Cl  (hounded  by  F ),  then  the  discretiza¬ 
tion  ofT  results  in  the  finite  dimensional  approximation  to  the  scattering  matrix,  to 
any  prescribed  precision. 


5.2  Recursive  Generation  of  Scattering  Matrices 

We  will  consider  a  case  when  a  compact  subset  D  of  domain  D  is  the  union  of  mutually 
disjoint  compact  sets  {A}  (see  Figure  5.3). 

It  turns  out  that  the  scattering  matrix  of  D  can  be  obtained  by  merging  the 
scattering  matrices  of  {D,}.  We  begin  with  introducing  the  requisite  notation.  Then 
we  present  the  merging  scheme  for  the  recursive  generation  of  scattering  matrices. 

5.2.1  Notation 

Suppose  that  A  =  {Fi,  12,  •  •  ■ ,  F^}  C  fl  is  a  set  of  closed  Jordan  curves.  Each  F,  G  A 
is  parameterized  by  its  length  7,  :  [0,  Z.,]  -+  R^,  and  Cli  C  Cl  '\s  the  region  bounded  by 
F,.  Suppose  further  that  for  1  <  i  <  m,  A  is  a  compact  subset  of  Cli,  function  G,: 
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Figure  5.4:  £),  C  fii  for  1  <  i  <  m 

fl,  X  n,  is  the  Green’s  function  for  domain  fi,,  and  function  iV,  :  [0,  L\  — >  R?  is 
the  interior  normal  to  F,. 

Assuming  the  domains  {Hi}  to  be  mutually  disjoint,  we  will  consider  the  compact 
subset  D  of  domain  17  defined  by  the  formula 

m 

c  =  u  ».■ 

i=l 

In  addition  to  operators  L,  P,  and  5  defined  by  (5.3),  (5.4),  and  (5.5)  in  the 
preceding  section,  we  will  require  the  operators  for  1  <  z  <  m 

I.. :  z:2(r.)  ^  L^{n,) 

P,  :  M{D,)  ^  lHr'^\d,) 
s,, :  M(D,)  ^  I^(r.) 

defined  by  formulae 


Lu{'^){x)  =  ^  /  ‘,^(7.(0) 

^  /I  Jo 


—  G.(x,7.(G)  ■  dt 


(5.9) 
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P.((t)(x)=  /  K[x,t)q{t)d<j{t)  iorx£R\D  (5.10) 

Jd, 

Si,{a){x)  =  f  K(x,t)q{t)da(t)  for  x  €  F,.  (5-11) 

JD, 

We  will  consider  equations  of  the  form 

P,{<r)  =  f  (5.12) 

with  /  e  L^{Di),  and  1  <  z  <  m. 

Definition  5.2  A  function  ifi  €  F^(r,)  will  be  called  total  incoming  potential  if 
for  any  x  G  A, 

P.(zTb.)(x)  =  L^ifQix),  (5.13) 

where  operators  Pi  and  La  are  defined  by  (5.10)  and  (5.9)  respectively,  and  a\D,  is 
the  restriction  of  the  charge  distribution  a  (defined  by  (5.4))  to  the  compact  subset 

A  c  D. 

Suppose  that  for  any  z  (1  <  z  <  m),  function  is  the  total  incoming  potential 
on  r,,  function  z/’,  is  the  outgoing  potential  induced  by  and  operator  a,  is  the 
scattering  matrix  for  the  domain  A.  Then 

a,  =  S„p-^L,„  (5.14) 

il),  =  ai-ifi,  (5.15) 

(see  (5.9),  (5.10),  (5.11),  (5.7),  and  (5.8)). 

We  will  also  require  operators 

L,  :  A(r)  ^  A(P), 

s, :  L’(r,)  ^  i'(r), 

Lj, :  L\r,) 
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defined  for  I  <  i,j  <  rn  via  formulae 

Li{ip){x)  =  —  dt  for  X  €  r„  (5.16) 

Sixhi  =  f  K{x,t)q{t)d<T{t)  for  X  6  r,  (5-17) 

Jd, 

Ljitf’i  =  f  K{x,t)q{t)d<r{t)  for  x  G  Fj.  (5.18) 

J  Di 

In  other  words,  the  operator  Si  converts  the  outgoing  potential  0,  on  F,  into  the 
scattering  potential  on  F,  and  the  operator  Lj,  converts  the  outgoing  potential  Wi  on 
Fi  into  the  scattering  potential  on  Fj. 


Definition  5.3  The  operai> 


S,  :  L\T) 


defined  by  the  formula 


1  I 

—  Li2Ct2  ••• 

-1 

(L,  \ 

—  L21Q1 

I 

L2 

~Lm2Ct2  ••• 

(5.19) 


will  be  referred  to  as  splitting  matrix,  provided  the  inverse  in  (5.19)  exists. 


5.2.2  Merging  Scheme  for  Scattering  Matrices 

Lemma  5.1  is  used  in  the  construction  of  the  merging  scheme  for  scattering  matrices, 
and  Theorem  5.2  gives  the  full  description  of  the  scheme. 
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Lemma  5.1  Suppose  that  is  an  incoming  potential  on  the  boundary  F,  and  <pi  is 
the  total  incoming  potential  on  F,  in  the  sense  of  (5.13).  If  the  inverse  in  (5.19) 
exists,  then 


<^2 


=  Sp-ip. 


(5.20) 


\V>m  / 


In  other  words,  the  splitting  matrix  maps  the  incoming  potential  p  on  T  to  the  total 
incoming  potentials  {tpi}  on  boundaries  {F^}. 


Proof:  For  any  J  (1  <  z  <  m),  the  total  incoming  potential  ipi  on  the  boundary  F, 
equals  to  the  sum  of  the  potential  from  F  and  the  scattering  potentials  from  {Fj} 
with  j  ^  i  and  1  <  i  <  m.  That  is, 

<Pi  =  Li<p-^Y^  Lijxfj.  (5.21) 

Combining  (5.15)  with  (5.21),  we  have 

ifi  =  Lig^  -f  Lijajifj.  (5.22) 

j#' 

Viewing  the  above  equations  as  a  m  x  m  linear  system,  we  obtain  (5.20).  ■ 


Theorem  5.2  (Recursive  Generation  of  Scattering  Matrices)  Given  scatter¬ 
ing  matrices  {oi}  for  domains  {Di},  the  scattering  matrix  a  of  domain  D  is  given  by 
the  formula 

a=^Siai  S2Q2  •••  SmCtm'^Sp,  (5.23) 

where  operators  {5i}  are  defined  by  (5.17),  and  the  splitting  matrix  Sp  is  defined  by 
(5.19). 


Proof:  Suppose  that  ip  is  an  incoming  potential  on  F,  and  cr  is  the  charge  distribution 
induced  by  p).  Then  by  definition  (see  equation  (5.1)),  the  induced  outgoing  potential 
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V’  assumes  the  form 


■ip{x)  =  f  K{x,t)  ■  q{t)  ■  da{t) 

Jd 

=  £  /  H{x,t)  ■  qit)  ■  dait) 


for  any  x  G  F. 

By  definition  of  the  operator  S',  in  (5.17), 

S,tl>,  —  f  K{x,t)q{t)dt, 
J  Di 

so  that  equation  (5.24)  assumes  the  form 


(5.24) 


(5.25) 


(5.26) 


Therefore, 


J]]  =  J] 

^  *^1^1  *5*20^2  '  *  '  ^ 

P2 

\  </^m  / 

Now,  the  conclusion  of  the  theorem  follows  from  the  combination  of  (5.27),  (5.7), 


and  Lemma  5.1.  ■ 
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Figure  5.5;  A  subset  D  of  Cantor  set  C“  in  domain  0 

5.3  Scattering  Matrices  in  Cantor  Sets 

Suppose  that  D  is  a  compact  subset  of  Cantor  set  C“,  and  the  set  D  is  enclosed  in  a 
domain  Q  with  its  boundary  denoted  by  F  (see  Figure  5.5).  The  boundary  F  will  be 
referred  to  as  frame  boundary. 

To  specifically  deal  with  the  problem  defined  in  section  2.2,  we  consider  scattering 
matrices  in  Cantor  sets,  with  Kix,t)  =  lnl|x  -  i||  and  q{t)  as  the  characteristic 
function  of  D.  Then  an  incoming  potential  is  harmonic  in  Cl,  and  an  outgoing 
potential  '5  is  harmonic  in  R^\D.  and  satisfying  the  far-field  condition  (2.8)  (see 
Remark  5.1). 

The  operators  P  :  L^{D)  —*  L^{R'^\D)  defined  by  formula  (5.4),  and  S  :  L^{D)  —>■ 
X^(F)  defined  by  formula  (5.5)  assume  the  form 


P(a){x)  =  f  In  ||x  —  t||c/cr(<), 

J  D 

(5.28) 

S{a){x)  =  [  \n\\x  —  t\\da it). 

JD 

(5.29) 
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In  this  section,  we  construct  an  analytical  apparatus  for  representing  scattering 
matrices  corresponding  to  subsets  of  Cantor  sets.  First,  we  discuss  the  representations 
of  incoming  and  outgoing  potentials  in  terms  of  single  and  double  layer  distributions. 

5.3.1  Representation  of  Potentials 

Lemma  5.2  below  is  an  obvious  statement  about  the  representations  of  incoming  and 
outgoing  potentials  in  terms  of  the  solutions  of  classical  boundary  value  problems  for 
the  Laplace  equation,  and  Lemma  5.3  is  a  statement  about  the  representations  of 
incoming  and  outgoing  potentials  in  terms  of  the  solutions  of  integral  equations  in 
the  potential  theory  (see  section  3.1).  Theorem  5.3  is  an  immediate  consequence  of 
Lemma  5.3.  and  Theorem  5.4  follows  immediately  from  Lemma  3.1  and  Theorem  5.3. 

Lemma  5.2  Suppose  that  .p  is  an  incoming  potential  on  the  boundary  F,  and  U  is 
an  outgoing  potential  on  F.  Then  the  incoming  and  outgoing  potentials  $  and  ^  are 
respectively  the  solutions  of  the  following  two  boundary  value  problems. 
fA.A)  Interior  Dirichlet  problem  (Incoming  Potential) 


A$(x)  =  0 
4>(x)  =  v?(7->(x)) 


for  X  E  Ll 
for  X  €  F 


(5.30) 


(BB)  Exterior  Neumann  problem  (Outgoing  Potential) 


A»l'(x)  =  0 
^'i>(x')  =  v(7“'(x)) 


for  X  ^  R^\Q 

for  X  G  F 


(5.31) 


with  'F  satisfying  the  far-field  condition  (2.8). 


Lemma  5.3  Suppose  that  'p>  is  an  incoming  potential  on  the  boundary  F,  and  xh  is  an 
outgoing  potential  on  F.  Suppose  further  that  a  dipole  distribution  pj  a.nd  a  charge 
distribution  p,  are  respectively  the  solutions  of  the  two  integral  equations, 


d 

rpji)  +  log  j|7(i)  -  7(i)||^j(i)d<  =  (^(: 


(5.32) 
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-  TTfisix)  +  log||7(x)-7(0||/i5(iM<  =  ^(x).  (5.33) 

Then  the  incoming  and  outgoing  potentials  can  he  represented  by  the  formulae 

^(2;)  =  ^  Q^^^og\\^{x)  - '){t)\\pait)dt,  (5.34) 

\I/{x)=  [  logl|7(i)-7(<)ll/i,(<)dt.  (5.35) 

Jo 

To  formulate  the  above  theorem  in  the  operator  notation,  we  define  four  operators 


Qd :  T"(r)  ^  I2(n), 

g.  :  L^T)  L\R\n), 

Pd :  I2(r)  I2(r), 

p, ;  L\T)  ^  L\r), 

via  formulae 

QdM{x)  =  j^  -^^logW-fix)  - -f{t)\\pd{t)dt,  (5.36) 

Q^{f^s)ix)=  I  log||7(x)-7(t)||//,(t)dt,  (5.37) 

Jo 

PdM{x)  =  Trpd{x)  + -Q^\og\\-f{x)-j{t)\\fid{t)dt,  (5.38) 

PsMix)  =  -Trpsix)  +  0^^  Jog  Il7(a^)  -  l/{l)\\fis[t)dt.  (5.39) 

Theorem  5.3  Suppose  that  (p  is  an  incoming  potential  on  the  boundary  F,  and  0 
is  an  outgoing  potential  on  F.  Then  the  incoming  and  outgoing  potentials  can  he 
represented  via  the  formulae 

$  =  g^p/V,  (5-40) 

^  =  Q.Pf^rh,  (5.41) 


where  operators  Qg,  Qd,  Ps  and  Pd  are  defined  by  (5.37),  (5.36),  (5.39),  and  (5.38) 
respectively. 
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Theorem  5.4  Suppose  that  L  :  L^(T)  L^(D)  is  the  operator  defined  by  the  formula 

(5.3).  Then 

L<p  =  iQdPf^‘p)\D,  (5.42) 

with  the  operators  Qd  and  Pd  defined  by  the  formulae  (5.36)  and  (5.38)  respectively. 

Remark  5.4  Under  certain  circumstances,  we  will  be  generating  scattering  matrices 
directly  by  using  their  definition  (see  (5.7)) 

a  =  SP-^L,  (5.43) 

where  operators  P,  S,  and  L  are  defined  by  formulae  (5.28),  (5.29),  and  (5-42)  re¬ 
spectively. 

5.3.2  Recursive  Generation  Of  Scattering  Matrices 

Suppose  that  {Z)i,D2,J93,jD4}  C  are  four  subsets  (boxes)  resulting  from  the 
subdivision  of  a  bigger  subset  (box),  and  the  set  D  is  the  union  the  four  subsets  {A} 

£>  =  U  A.  (5.44) 

«=i 

Suppose  further  that  D  is  enclosed  in  s  square  O  with  its  boundary  denoted  by  F. 
The  square  SI  will  be  referred  to  as  frame  domain  (box)  while  F  is  referred  to  as  frame 
boundary. 

Suppose  that  for  any  integer  i  (1  <  i  <  m),  the  set  A  is  enclosed  in  a  square  fl, 
with  its  boundary  denoted  by  F,  (see  Figure  5.6).  Within  the  tree  structure  of  the 
Cantor  set  C“  (see  section  2.1),  we  will  refer  the  frame  boxes  of  neighbor  boxes  as 
frame  neighbor  boxes,  and  the  frame  box  of  a  parent  box  as  parent  frame  box. 

In  this  section,  we  obtain  the  scattering  matrix  a  for  D  from  scattering  matrices 
{Q;i,a2,Qr3,Q;4}  for  domains  {I?i,A?  A,  A}-  First,  we  need  the  representations  of 
operators 


L,  :  A(F)  ^  A(F.), 
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Figure  5.6:  Four  Subsets  of  C“  and  Their  Frame  Boxes 

5,  :  L\Ti)  L\T), 

L,,  :  L\r.)  -  L^{Tj), 

defined  by  formulcie  (5.16),  (5.17),  and  (5.18)  respectively. 

Theorems  5.5  and  5.6  below  follow  from  Theorem  5.3  immediately.  Lemma  5.4 
is  the  immediate  consequence  of  Lemma  5.1,  and  Theorem  5.7  is  the  consequence  of 
Theorem  5.2  and  Lemma  5.4.  Theorems  5.5  and  5.6  describe  representations  of  the 
operators  5,,  and  defined  by  formulae  (5.16),  (5.17),  and  (5.18).  And  Theorem 
5.7  describes  the  merging  scheme  for  the  recursive  generation  of  scattering  matrices 
for  subsets  of  Cantor  sets. 

Theorem  5.5  Suppose  that  Li  :  L^(r)  L^(rj)  is  the  operator  defined  via  formula 
(5.16).  Then 

Liip  =  (Q<iF, -V)lr,  (5.45) 


with  operators  Qi  and  Pd  defined  by  formulae  (5.36)and  (5.38)  respectively. 
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Similar  to  operators  Qs  defined  by  formula  (5.37)  and  P*  defined  by  formula  (5.39) 
for  domain  fl,  we  define  operators  on  domain  fl,  for  1  <  i  <  4, 

g,,  :  L\T,)  L\R'^\Q,), 

P,,  :  L\r,)  ^  L^(r.), 

by  formulae 

QsA^s){x)  ^  f  logll7.(a:) -7i(t)||(Ts(t)dt,  (5.46) 

JQ 

Ps,i{<T^){x)  =  -7r<j,(i)  +  log  il7i(2:)  -  7,(t)llcr,(t)<it.  (5.47) 

Theorem  5.6  Suppose  that  operator 

Si :  Z^(r.)  x^(r), 

is  defined  via  formula  (5.17),  and  operator 

Ljr.  L\Ti)  ^  L^Tj)  fori:^j, 
is  defined  via  formula  (5.18).  Then 


(5.48) 

Ljii’  =  {Qs,iPfl^)\rj, 

(5.49) 

with  operators  Q,  i  and  P,^i  defined  by  formulae  (5-46)  and  (5.47)  respectively. 
Remark  5.4  The  splittering  matrix 

'  i"(r,) ' 


s, :  i^(r) 
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is  given  by  the  expression  (see  formula  (5.19)) 


I  —L12012  — L13O3  —L\ia^  ^ 

-1 

(lA 

—  ^2101  I  — L23O3  — X24O4 

£2 

—  Lz\OCl  —£320:2  I  ~I'2iOC4 

£3 

\  —£4101  —L420C2  —£4303  I  j 

<  -^4  j 

(5.50) 


with  Li  defined  by  (5-45)  for  1  <  i  <  4,  and  Lij  defined  by  (5-49)  for  1  <  i,j  <  4. 

Lemma  5.4  Suppose  that  v?  is  an  incoming  potential  on  the  boundary  F,  and  for 
1  <  I  <  4,  is  the  total  incoming  potential  on  F,-  in  the  sense  of  (5.13).  Then 


'  ¥>1  '' 
(^3 


(5.51) 


with  the  splitting  matrix  Sp  defined  by  formula  (5.50). 

In  other  words,  the  splitting  matrix  maps  the  incoming  potential  (p  on  T  to  the 
total  incoming  potentials  {<^,}  on  boundaries  {F,}. 


Theorem  5.7  Suppose  that  for  1  <  i  <  4,  the  scattering  matrix  for  the  compact  set 
Di  is  denoted  by  o;.  Then  the  scattering  matrix  a  for  the  set  D  =  Uf_i  Di  is  given 
by  the  merging  formula 

^  Sioii  5202  Szctz  54Q4  ^  5p,  (5.52) 

with  operator  Si  defined  by  formula  (5-48)  for  1  <  z  <  4,  and  the  splitting  matrix  Sp 
defined  by  formula  (5.50). 


Remark  5.5  Due  to  the  self-similarity  in  Cantor  sets  (see  section  2.1),  we  only  need 
to  compute  one  scattering  matrix  per  level. 
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Figure  5.7:  Frame  Boxes  for  Cantor  Set  C° 

5.3.3  Discretization  of  Scattering  Matrices 

From  the  preceding  sections,  it  is  clear  that  the  construction  of  scattering  matrices  for 
subsets  of  Cantor  sets,  either  directly  or  recursively,  depends  on  the  choice  of  frame 
boundaries  (or  frame  boxes)  (see  Remark  5.4  and  Theorem  5.7).  For  a  Cantor  set 
C“,  we  recursively  generate  frame  boxes  for  subsets  of  C“  such  that  frame  boxes  at 
the  same  level  are  mutually  disjoint,  and  the  distance  between  a  box  and  its  frame 
box  equads  to  the  distance  between  two  neighbor  frame  boxes  (see  Figure  5.7). 

With  the  above  choice  of  frame  boxes,  we  first  represent  incoming  and  outgoing 
potentials  $  and  ^  numerically  in  the  forms  described  in  Theorem  5.3,  in  order  to 
discretize  scattering  matrices  for  subsets  of  Cantor  sets. 

Lemma  5.5  below  introduces  a  set  of  orthogonal  polynomials,  and  Lemmas  5.6 
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eind  5.59  describe  the  two-dimensional  integrations  involved  in  the  application  of  the 
Galerkin  method  to  the  numerical  representation  of  incoming  and  outgoing  potentials 
(see  sections  3.2  and  5.3.1). 


Lemma  5.5  Suppose  that  {Pi{x),  P2ix),- ■  ■ ,  Pn{x),- •  ■}  is  a  set  of  polynomials  de¬ 
fined  by  the  formulae 


(5.53) 


Then 


J  Pi{x)Pj{x)dx  =  0  for  i  ^  j. 


(5.54) 


The  orthogonal  polynomials  defined  above  are  the  well-known  Legendre  poly¬ 
nomials.  The  roots  of  Pn{x)  will  be  referred  to  as  the  Legendre  nodes.  We  will 
use  the  Legendre  polynomials  in  the  Galerkin  method  as  the  basis  functions  to  rep¬ 
resent  incoming  and  outgoing  potentials  numerically  in  terms  of  the  representations 
described  in  Theorem  5.3. 


Lemma  5.6  Suppose  that  Q,  is  a  square  of  size  two  with  its  boundary  denoted  by  F, 
and  {Pi(x),  P2{x),  •  •  • ,  Pn{x),  *  •  •}  are  the  Legendre  orthogonal  polynomials  defined  in 
Lemma  5.5.  Suppose  further  that  p  is  either  the  solution  of  equation  (5.32)  or  of 
equation  (5.33),  and  the  solution  p  is  approximated  by  the  truncated  expansion  on 
each  side  of  the  square  ft: 

=  (5.55) 

i=0 

for  each  of  the  two  horizontal  sides  of  U,  and 

Mn(y)  =  ^T;iP,(y),  (5.56) 

«=o 

for  each  of  the  two  vertical  sides  of  fi.  Then  there  are  only  two  types  of  integrations 
involved  in  the  matrices  of  the  type  defined  by  formula  (3.17)  in  the  Galerkin  method 
described  in  section  3.2.  The  two  integrations  are  of  the  form 
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for  any  two  adjacent  sides  of  domain  Q,,  and 

for  non-adjacent  sides  of  domain  Ll,  where  1  <  i,j  <  n. 

The  following  lemma  is  obtained  immediately  by  writing  down  the  integral  /,j 
(defined  by  (5.57))  in  polar  coordinates. 

Lemma  5.7  Suppose  that  for  1  <  i,j  <  n,  integration  Lj  is  defined  by  (5.57).  Then 
r~  r-  r— 

/,j  =  { y  ^  y  A  J  ^  y  ’  }P,(r  cos^  —  l)Pj(rsin^  —  l)sin^dr<i0.  (5.59) 

Remark  5.6  Integrations  7,^  defined  by  formula  (5.59)  and  Jij  defined  by  formula 
(5.58)  for  1  <  2,  J  <  n,  are  now  integrals  of  smooth  functions.  They  can  be  computed 
by  the  Gaussian  quadrature  rule  based  on  Legendre  nodes  (see,  for  example,  [50]). 
Thus,  operators  Pd  and  P,  can  be  represented  numerically  by  the  formulae  (5.38)  and 
(5.39)  respectively. 

On  the  other  hand,  with  the  approximated  solution  pn  (defined  by  (5.55)  and 
(5.56))  to  p  the  solution  of  either  equation  (5.32)  or  (5.33),  operators  Qd  and  Qs 
defined  by  formulae  (5.36)  and  (5.37)  can  be  represented  numerically  by  using  the 
Gaussian  quadrature  rule. 

Therefore,  incoming  and  outgoing  potentials  $  and  ^  are  represented  numerically 
by  formulae  (5.40)  and  (5.41)  respectively. 

Beised  on  the  above  choice  of  frame  boxes,  Tables  5.1  and  5.2  list  the  number  of 
Legendre  nodes  needed  on  frame  boundaries  for  the  representation  of  incoming  and 
outgoing  potentials  to  single  and  double  precision  respectively. 
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e  =  10  ^  (Absolute  Error) 


Ratio  a 

0.1 

0.2 

0.3 

0.35 

0.4 

0.45 

Number  of  Nodes 
Per  Side 

12 

18 

30 

46 

66 

100 

Table  5.1:  Number  of  Legendre  Nodes  for  Single  Precision 


c  =  10  (Absolute  Error) 


Ratio  a 

0.1 

0.2 

0.3 

0.35 

0.4 

0.45 

Number  of  Nodes 
Per  Side 

30 

44 

60 

80 

120 

200 

Table  5.2:  Number  of  Legendre  Nodes  for  Double  Precision 

Remark  5.7  Two  observations  can  be  easily  made  from  the  above  two  tables. 

1.  The  number  of  Legendre  nodes  needed  to  obtain  double  precision  is  only  twice 
the  number  of  nodes  needed  for  single  precision. 

2.  The  number  of  nodes  needed  increases  rapidly  with  the  increase  in  ratio  a. 

With  the  frame  boundaries  described  above,  we  are  ready  to  discretize  scattering 
matrices  for  subsets  of  Cantor  sets. 

Suppose  the  the  set  of  points  2:27  •  •  • ,  is  the  approximation  to  the  compact 
set  D  C  and  the  set  of  points  {xi,a:2,- •  •  ,Xp}  are  the  Legendre  nodes  on  the 
boundary  P. 

The  operator  P  defined  by  formula  (5.28)  is  discretized  in  a  manner  similar  to 
that  of  the  integral  equation  (2.9).  In  other  words,  the  discretization  of  P  is  a  matrix 
P  defined  by  the  formulate 

{P),j  =  In \\zi  -  Zj\\  for  i  j  and  1  <  i,j  <  m,  (5.60) 


{P)ii  =  I  ■  logo  +  ^  for  1  <  i  <  m, 


(5.61) 
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where  6  is  a  contant  (see  (2.11)). 

Similarly,  the  operator  S  defined  by  (5.29)  is  discretized  as  a  matrix  5  defined  by 
the  formulae 


(S),j  =  In  ||x,-  —  ZjW  for  1  <  2  <  p  and  1  <  j  <  m.  (5.62) 

The  operator  Qd  defined  by  (5.36)  and  the  operator  defined  by  (5.46),  are 
discretized  by  the  Gaussian  quadrature  rule  based  on  Legendre  nodes.  The  operator 
Pd  defined  by  (5.38)  and  the  operator  defined  by  (5.47)  are  discretized  by  the 
Galerkin  method  described  in  section  3.2. 

Remark  5.8  Given  a  discretization  of  operators  P,  S,  Pd,  Qd,  Qa,i  o.nd  Ps,i  (defined 
by  the  formulae  (5.28),  (5.29),  (5.38),  (5.36),  (5.46)  and  (5.41)  respectively) ,  all  the 
other  operators  L,  L,,  Si,  Lij  and  Sp  for  1  <  i,j  <  4,  are  represented  numerically 
by  the  formulae  (5.42),  (5.45)  ,  (5.48),  (5.49)  and  (5.50).  Thus,  scattering  matrices 
can  be  computed  either  directly  by  using  formula  (5.43),  or  recursively  by  using  the 
merging  scheme  described  in  Theorem  5.2. 


Chapter  6 


The  Fast  Direct  Algorithm 


In  this  chapter,  we  describe  a  direct  algorithm  for  the  rapid  solution  of  the  Laplace 
equation  on  regions  with  fractal  boundaries.  The  algorithm  exploits  the  fact  that 
for  any  given  ratio  a,  interactions  at  any  level  in  the  Cantor  set  C“  are  of  low  rank 
(the  ranks  depend  only  the  constant  ratio  a  for  generating  the  Cantor  set,  and  do 
not  depend  on  the  sizes  of  boxes  and  number  of  points  inside).  The  low  rank  of 
interactions  is  reflected  in  the  coelBcient  matrix  in  equation  (2.10)  as  the  low-rank 
of  its  oflf-diagonal  submatrices  (see  Figure  6.1).  Thus,  we  can  recursively  compress 
these  matrices  of  low  rank  without  actually  generating  them. 

To  be  more  specific,  let  us  consider  four  subsets  (boxes)  in  a  Cantor  set  C“, 
depicted  in  Figure  6.2.  They  are  boxes  of  size  d.  and  the  distance  between  any  two  of 
them  is  (1  —  2a)d.  The  interactions  between  them  are  of  low  rank  (see  section  4.2), 
and  can  be  represented  via  scattering  matrices  (see  section  5.3). 

Starting  with  the  hierarchical  structure  of  a  Cantor  set  C“  (see  section  2.1),  we 
proceed  by  introducing  a  set  of  frame  boxes  arranged  in  a  tree  structure  (see  section 
5.3).  For  a  given  precision  e,  we  determine  the  number  of  Legendre  nodes  needed  on 
frame  boundaries  for  the  representation  of  potentials  (see  Tables  5.1  and  5.2).  Then 
we  precompute  the  inverses  of  operators  Pd  and  defined  by  formulae  (5.38)  and 
(5.39),  via  the  classical  Galerkin  method  (see  section  5.3  1). 
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Figure  6.1:  Approximation  of  Cantor  set  and  Coefficient  Matrix 

To  describe  the  algorithm,  we  need  the  following  notation. 

L  Number  of  levels  in  the  approximation  of  the  Cantor  set  C“. 

N  N  =  the  size  of  the  approximation  of  the  Cantor  set  C“ 

p  The  number  of  Legendre  nodes  on  each  side  of  frame  boundaries 

for  the  representation  of  potentials  to  a  given  precision  e 
Li  The  number  of  a  level  on  which  a  scattering  matrix  is  computed  directly. 

JTi  =  4^“^^ ,  size  of  linear  systems  to  be  solved  directly 

Af  The  m  X  m  matrix  Aj  is  the  restriction  of  matrix  A  in  (2.11)  onto  a 
subset  {ibox)  at  level  Li.  In  other  words,  (A/),j  =  In  ||z,  —  Zj||,  and 
{Aj)ii  =  Llna  +  (5,  where  points  {zj,  •  •  • ,  Zn]  C  ibox,  and  (5  is  a  constant. 
The  fast  direct  algorithm  is  a  two-pass  procedure.  In  the  first  (bottom-up)  pass, 
we  compute  the  scattering  matrix  for  level  Li  directly  by  using  formula  (5.4.3),  and 
scattering  matrices  for  all  coarser  levels  (level  number  <  )  by  using  the  merging 

scheme  described  in  Theorem  5.7.  In  the  second  (top-down)  pass,  we  generate  total 
incoming  potentials  on  frame  boundaries  up  to  level  Li  by  using  Lemma  5.4.  Finally, 
we  solve  4^'  small-scale  linear  systems  of  size  m  x  m  directly  at  level  L\ . 

Following  is  a  formal  description  of  the  algorithm. 
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(1  -2a)j| 


Figure  6.2:  Four  Subsets  of  Cantor  set  C“ 

The  Algorithm 

Initializaiton 

Comment  [Computations  in  the  initialization  are  done  once  for  all.] 

Step  1 

Comment  [  Given  a  real  number  a  (0  <  a  <  |)  and  an  integer  L,  construct  the 
approximation  of  Cantor  set  C“,  and  its  frame  boxes] 
do  lev  =  0.1,2,  -  ■■  ,L 

do  ibox  =  1, 2,  •  •  • , 

Divide  each  box  into  four  corner  boxes  according  to  the  constant  a. 
Construct  the  frame  box  for  ibox  (section  5.3). 

endo 

endo 

do  lev  =  L 

do  ibox  ~  1, 2.  •  ■  • ,  4^ 

Compute  the  center  of  ibox. 
enddo 


enddo 
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Step  2 

Comnienti  For  a  given  precision  t.  precompute  operators  .  Pp^  •  and  .4^’] 
do 

Determine  number  of  Legendre  nodes  p  (see  Tables  5.1  and  5.2) 

Generate  the  inverses  of  operators  Pd  and  P,  defined  in  (5.38)  and  (5.39),  via 
the  classical  Clalerkin  method. 

Compute  the  inverse  of  m  x  rn  matrix  Aj  directly. 

enddo 

Upward  Pass 

CommentiCompute  scattering  matrices  and  splitting  matrices] 

Step  3 

do  lev  =  Li 

Compute  the  discretized  operators  L  and  5  defined  by  (5.42)  and  (5.29). 
Compute  the  scattering  matrix  directly  via  formula  (5.43):  a  =  SAJ^L. 

endo 

Step  4 

do  lev  =  L\  —  \  .  L\  -  2.  •  •  •  .  1. 0 
do  z  =  1,  •  ■  ■ .  4 

Compute  operators  T,  and  5,  (defined  by  formulae  (5.45)  and  (5.48)) 
via  the  Gaussian  quadrature  rule  based  on  Legendre  nodes, 
do  j  =  1.  ■  ■  • ,  4 

Compute  Lij  defined  by  (5.49)  via  the  Gaussian  quadrature  rule. 

endo 

endo 

Compute  the  splitting  matrix  5p  by  using  formula  (5.50). 

Compute  the  scattering  matrix  a  via  the  merging  scheme  in  Theorem  5.7. 


endo 
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Downward  Pass 

Comment  [Splitting  matrices  are  now  available.  Compute  total  incoming  potentials 
on  all  frame  boundaries  up  to  level  Li] 


Step  5 


do  lev  =  1,2,  -  -  •  ,1.1 

do  ibox  =  1, 2,  •  •  ■ ,  4^®" 

Compute  total  incoming  potential  by  means  of  Lemma  5.4. 
enddo 

endo 


Step  6 


do  ibox  =  1,2,  •  •  •  ,4^1 

Solve  m  X  m  linear  system  directly  by  computing 
where  operator  L  is  computed  at  Step  3. 

endo 

Remark  6.1  Suppose  that  ibox  is  a  fixed  box  at  level  1.  Then  in  the  splitting  pro¬ 
cess,  the  total  incoming  potential  on  the  the  frame  boundary  of  ibox  can  be  computed 
independently  from  those  on  the  other  frame  boundaries  of  boxes  at  the  same  level. 

Thus,  we  can  obtain  a  part  of  the  solution  independently  from  the  rest  of  the 
solution  if  only  a  part  of  the  solution  is  desired. 


A  brief  analysis  of  the  algorithmic  complexity  is  given  below. 
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Step  Number 

Operation  Count 

Explanation 

1 

OiN) 

AN  boxes  (squares)  are  involved.  Each 
box  is  determined  by  its  center  and  size. 

2 

0(m3  +  p3) 

Operators  Pd  and  P^  are  of  size  4p  x  4p, 
and  the  Operator  ^4/  is  of  size  m  x  m. 

3 

0{mp  +  m^p) 

Operator  L  is  of  size  m  x  4p. 

Operator  S  is  of  size  Ap  x  m. 

Operator  A^^  is  of  size  m  x  m. 

4 

0{p^  log  N) 

Operator  5,,  and  Lij  is  of  size  4p  x  4p. 

Sp  is  of  size  16p  x  4p.  Operator  a  is  of 
size  4p  X  4p.  There  are  log  N  levels. 

5 

0{p^N) 

The  computation  of  the  total  incoming 

potential  on  each  frame  boundary 
requires  p^  operations.  There  are 
(<  N)  frame  boundaries  involved. 

6 

0{m^pN) 

Operator  A~j^  is  of  size  m  x  m.  Operator 

L  is  of  size  m  x  4p.  Potential  (p,  is  a 

vector  of  size  4p.  Computations  of 
are  done  4^’(<  N)  times. 
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The  time  complexity  of  the  algorithm  is  therefore 

(/3iP^  +  /32m^p)  •  N  +  03p^  •  log  N,  (6.1) 

where  the  constant  m  is  normally  chosen  to  be  m  ~  256,  and  the  constant  p  depends 
on  the  geometry  of  a  given  fractal  boundary,  and  the  choice  of  frame  boundaries  (see 
Tables  5.1  and  5.2),  and  the  constants  /?2,  and  /?3  depend  on  the  computer  system, 
implementation,  lajiguage,  etc. 

Remark  6.2  Given  scattering  matrices  and  total  incoming  potentials,  the  evaluation 
of  the  potential 

^(ar)  =  /  In ||x  —  t||d(7(<)  (6.2) 

at  any  point  x  €  R?\C°’  requires  at  most  O(log  N)  operations,  where  a  is  the  charge 
distribution  over  C“. 


Chapter  7 


Numerical  Experiments 


We  have  implemented  the  fast  direct  algorithm  of  the  preceding  section  in  Fortran 
77.  The  program  is  capable  of  computing  either  whole  or  part  of  the  solution,  and  of 
evaluating  potential  at  any  point.  We  used  our  algorithm  to  compute  the  harmonic 
measure  on  Cantor  sets  with  the  dimension  of  support  of  harmonic  measure  deter¬ 
mined  via  the  entropy  of  a  set  of  charges  distributed  on  a  Cantor  set,  the  charges 
representing  the  solution  of  the  Laplace  equation  with  unity  as  the  boundary  condi¬ 
tion  on  that  Cantor  set. 

We  will  need  the  following  terminology  to  describe  our  numerical  experiments. 

1.  HausdorfF  dimension  of  Cantor  set  C“  is  given  by  the  formula 

=  log4/log(-).  (7.1) 

a 

2.  Dimension  of  Support  of  harmonic  meaisure  on  C“  is  given  by  the  formula 

d  =  /f/log(-),  (7.2) 

a 

where  H  is  the  entropy  of  the  system. 


3.  Approximation  of  entropy  H 
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where 


1  ^ 

fff,  =  ^(7,  -loglTj 


(7.4) 


1=1 


with  {(T,}  the  scaled  charge  distribution  over  the  system  ==  !)• 


We  considered  fractal  boundaries  of  Cantor  type  with  ratio  0.1,  0.3,  and  0.45. 
For  the  first  two  experiments  (a  =  0.1, 0.3),  the  potentials  on  frame  boundaries  are 
represented  to  double  precision  while  they  axe  represented  to  10  digits  for  the  third 
experiment  (a  =  0.45).  The  size  of  linear  systems  inverted  directly  at  the  final  stage 
ha.s  been  chosen  to  be  m  =  256.  The  constant  6  in  the  coeflBcient  matrix  A  defined 
by  formula  (2.11)  is  chosen  to  be 


S/q  +  {y-  yoy)dxdy 

fSadidy 


(7.5) 


where  Q  is  a  square  centered  at  (xo,yo)  on  the  finest  level  of  the  recursive  generation 
of  L-th  approximation  of  Cantor  set  C“  (see  section  2.1).  All  calculations  have  been 
conducted  on  a  Sparc  II  workstation. 

The  results  are  summarized  in  tables  7.1,  7.2,  7.3,  7.2,  7.5,  and  7.6.  In  tables  7.1, 
7.3,  and  7.5,  the  first  column  is  the  size  of  the  approximation  to  a  Cantor  set.  The 
second  column  is  the  number  of  levels  in  the  generation  of  a  Cantor  set.  The  third 
column  is  the  actual  CPU  time  of  the  fast  direct  algorithm  of  the  preceding  section. 
The  forth  column  is  either  the  CPU  time  or  estimated  CPU  time  of  the  combined 
algorithm:  the  Conjugate  Gradient  (CG)  algorithm  combined  with  the  Fast  Multipole 
Method  (FMM)  (see  [46]).  The  leist  column  is  the  estimated  timing  for  the  Gaussian 
Elimination  (of  course,  it  is  given  here  only  for  comparison  purposes). 

The  following  observations  can  be  made  from  Tables  7.1,  7.3,  and  7.5. 


1.  Although  our  fast  algorithm  asymptotically  requires  0{N)  operations,  the  ac¬ 
tual  running  time  of  the  algorithm  as  observed  from  the  numerical  experiments 
seems  to  behave  like  log  N,  due  to  the  fact  that  the  constant  ^3  in  formula  (6.1) 
is  rather  large  compared  to  the  constants  and  02-  The  constant  0^  in  formula 
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(6.1)  can  be  substantially  reduced  either  by  a  better  choice  of  frame  boundaries, 
or  by  improving  the  representation  of  potentials  . 

2.  For  any  N  >  4096,  our  algorithm  is  faster  than  the  combined  algorithm  (the 
CG  method  combined  with  the  FMM,  see  [46]). 


3.  The  performance  of  our  algorithm  deteriorates  with  the  increase  in  ratio  a  as  is 
expected  (see  Tables  5.1  and  5.2). 


Tables  7.2,  7.4,  and  7.6  summarize  some  of  our  numerical  experiments  in  which 
the  dimension  of  support  of  harmonic  measure  was  computed.  The  part  A  of  the 
tables  is  the  computed  dimension  of  the  support  d.  The  part  B  of  the  tables  is  the 
the  Richarson  extrapolation  of  the  first  order  in  terms  of  levels  by  the  formula 

^  Tldji 

m  —  n 


(7.6) 


with  dm  and  d„  the  dimensions  of  support  at  levels  m  and  n  respectively.  The  part 
C  of  the  tables  is  the  Richarson  extrapolation  of  the  second  order  in  terms  of  levels 
defined  by  the  formula 


d  = 


rrrdm,i  -  rrdim 


(7.7) 


m‘’  — 

with  dm,i  and  the  extrapolated  dimensions  of  support  defined  by  formula  (7.6). 
The  following  observations  can  be  made  from  Tables  7.2,  7.4,  and  7.6. 


1.  The  dimension  of  the  support  of  harmonic  measure  converges  linearly  with  the 
number  of  levels. 


2.  Despite  numerical  problems  often  associated  with  the  type  of  Richarson  ex¬ 
trapolation  described  above  (see  [50]),  the  experiments  did  show  rea.sonable 
convergence  rates:  the  first  order  extrapolation  in  Table  7.4  is  consistent  with 
Carleson  theorem  (see  [12]).  However,  the  situation  has  to  be  analyzed  carefully. 

3.  Even  though  our  fast  algorithm  permits  simulations  on  a  very  large  scale  to  be 
performed,  the  field  being  investigated  clearly  could  benefit  from  the  application 
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of  supercomputers,  and  possibly  parallel  machines.  Investigation  of  parallel 
implementations  of  algorithms  of  the  type  described  above  will  be  reported  at 
a  later  date. 

For  illustration,  the  dimension  of  support  is  plotted  as  the  function  of  the  Haus- 
dorff  dimension  in  Figures  7.4,  7.5,  7.7,  and  7.7.  Charge  distributions  are  also  plotted 
in  Figures  7.8,  7.9,  and  7.10  along  the  horizontal  lines  on  which  the  charges  are 
located. 

Remark  7.1  The  value  (7.5)  of  the  constant  6  has  been  chosen  empirically,  and  no 
claim  is  made  here  as  to  its  optimality. 

The  following  is  a  recapitulation  of  the  other  notation  to  be  used  in  the  illustration 
of  our  numerical  experiments. 

a  —  ratio  for  generating  a  Cantor  set. 

p  —  number  of  Legendre  nodes  on  each  side  of  frame  boundaries. 
m  —  size  of  linear  systems  inverted  directly. 

e  —  precision  to  which  incoming  and  outgoing  potentials  are  represented. 

L  —  number  of  levels  in  the  approximation  of  a  Cantor  set. 

—  size  of  the  linear  system  (2.10)  to  be  solved,  N  =  4^. 
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Figure  7.1:  Cantor  Set  with  a  =  0.1 

Example  7.1  Harmonic  Measure  on  Cantor  Set  C“  with  Ratio  a  =  0.1. 


a  =  0.1,  p  =  30,  m  =  256,  e  =  10 


N 

Levels 

ro/g(  Minutes) 

Tcg&fm  Af  (Hours) 

Tge  (Estimated) 

4,096 

6 

6 

0.4 

19.1  Hours 

16,384 

7 

9 

3.3 

51  Days 

65,536 

8 

11 

26.4  (est.) 

9  years 

262,144 

9 

14 

211.2  (est.) 

572  years 

1,048,576 

10 

19 

1689.6  (est.) 

366283  years 

Table  7.1;  Comparison  of  Timings  (a  =  0.1) 
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Hausdorff  Dimension  D  =  0.602059991 327962£'  +  00 


Level  Number  L 

Dimension  of  Support  d 

6 

0.599329045439459E+00 

7 

0.599189754023720E+00 

8 

0.599085067779161E+00 

9 

0.599003595888791E+00 

10 

0.5989384071 16325E+00 

First  Order  Extrapolation 


Level  Number  L 

Dimension  of  Support  d 

6,7 

0.598354005529286E+06 

7,8 

0.598352264067248E+00 

8,9 

0.598351820765831E+00 

9,10 

0.598351708164130E+00 

Second  Order  Extrapolation 


Level  Number  L 

Dimension  of  Support  d 

6,7,8 

0.59834703968 1135E+00 

7,8,9 

0.59835026921 0872E+00 

8,9,10 

0.598351257757325E+00 

Table  7.2;  Dimension  of  Support  d  on  Cantor  Set  C“  with  a  =  0.1 
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Figure  7.2:  Cantor  Set  with  a  =  0.3 

Example  7.2  Harmonic  Measure  on  Cantor  Set  C“  with  Ratio  a  —  0.3. 


a  =  0.3,  p  =  60,  m  =  256,  c  =  10 


Table  7.3:  Comparison  of  Timings  (a  =  0.3) 
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Hausdorff  Dimension  D  =  0.11ol43328498689£'  +  01 


Level  Number  L 

Dimension  of  Support  d 

6 

0.102592656926789E+01 

7 

0.101715537867085E+01 

8 

0.101051673123065E+01 

9 

0.100533476169737E+01 

10 

0.100118336944343E+01 

First  Order  Extrapolation 


Level  Number  L 

Dimension  of  Support  d 

6,7 

0.964528235088610E+00 

7,8 

0.964046199149250E+00 

8,9 

0.963879005431 130E+00 

9,10 

0.963820839157970E+00 

Second  Order  Extrapolation 


Level  Number  L 

Dimension  of  Support  d 

6,7,8 

0.962600091331 170E+00 

7,8,9 

0.963293827417710E-I-00 

8,9,10 

0.963588174065331E+00 

Table  7.4:  Dimension  of  Support  d  on  Cantor  Set  C“  with  a  =  0.3 
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Figure  7.3:  Cantor  Set  with  a  =  0.45 

Example  7.3  Harmonic  Measure  on  Cantor  Set  with  Ratio  a  =  0.45. 


a  =  0.45,  p  =  80,  m  =  256,  e  =  10 


N 

Levels 

Ta/5(  Minutes) 

TcGiiFMM{HoUTs) 

Toe  (Estimated) 

4,096 

6 

122 

1.6 

19.1  Hours 

16,384 

7 

181 

13.0 

51  Days 

65,536 

8 

240 

103.7(est.) 

9  years 

262,144 

9 

297 

829.3(est.) 

572  years 

1,048,576 

10 

360 

6634.2(est) 

366283  years 

Table  7.5:  Comparison  of  Timings  (o  =  0.45) 
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HausdorfF  Dimension  D  =  0.173610644917543£  +  01 


Level  Numbers 

Dimension  of  Support 

6 

0.122269892939995E+01 

7 

0.119171383040669E+01 

8 

0.116808624449636E+01 

9 

0.1 14954943391 115E+01 

10 

0.113465307668069E+01 

First  Order  Extrapolation 


Level  Number  L 

Dimension  of  Support  d 

6,7 

0.100580323644713E-f01 

7,8 

0.100269314312405E+01 

8,9 

0.100125494922947E+01 

9,10 

0.100058586160655E+01 

Second  Order  Extrapolation 


Level  Number  L 

Dimension  of  Support  d 

6,7,8 

0.9933628631 54809E+00 

7,8,9 

0.996221270598445E+00 

8,9,10 

0.99790951 1114859E+00 

Table  7.6:  Dimension  of  Support  d  on  Cantor  Set  C“  with  a  =  0.45 


X 


/! 

3.50  1.00  1.50 


7.4;  Dimension  of  Support  d  vs.  Hausdorff  Dimension  D 
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Nlev=6.7 

Niev=7.’8 

Klev=8.9 

Mev=^.10 


X 


Figiire  7.5:  Dimension  of  Support  d  vs.  Hausdorff  Dimension  D 
(First  Order  Extrapolation) 
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Figure  7.6:  Dimension  of  Support  d  vs.  HausdorfF  Dimension  I) 
(Second  Order  Extrapolation) 
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Y 


Nlev=6 

Nie7=7 . 

Nlev=8 

Nlev=10 

Nlev=6.7“  ■ 

Nle^73 

Niev“8.^ 

Nlev=9.10 

Niev=8.9.r0 

Nlev=7?8‘.9'' 

Nrev=67778’ 


X 


Figure  7.7:  Dimension  of  Support  d  vs.  HausdorfF  Dimension  D 
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Figure  7.8:  Charge  Distribution  for  a  =  0.1  {N  =  4096) 


Chapter  8 


Conclusions 


and  Generalizations 


We  have  presented  an  0{N)  direct  algorithm  for  the  rapid  solution  of  the  Laplace 
equation  on  regions  with  fractal  boundaries.  In  the  algorithm,  operators  of  low  rank 
are  recursively  compressed,  and  the  inverse  is  constructed  in  a  compressed  form  so 
that  it  can  be  applied  to  a  vector  rapidly.  The  algorithm  is  capable  of  generating  only 
a  part  of  the  solution  if  desired.  The  evaluation  of  potential  at  any  point  requires 
O(logfV)  operations. 

The  fast  direct  algorithm  of  this  thesis  admits  far-reaching  generadizations.  Fol¬ 
lowing  are  some  of  the  examples. 

Harmonic  Measure  On  Cantor  Sets  In  Three  Dimensions 

It  is  straightforward  to  generalize  the  algorithm  of  this  thesis  to  solve  the 
Laplace  equation  in  three  dimensions  on  regions  with  fractal  boundaries.  In 
R^,  the  behavior  of  harmonic  measure  for  Cantor  sets  is  completely  unknown. 
Peter  Jones  recently  raised  a  question  about  the  determination  of  the  actual  val¬ 
ues  of  the  dimension  of  the  support  of  harmonic  measure  on  fractals  of  certain 
types  in  and  conjectured  that  the  dimension  of  the  support  in  should  be 
always  less  than  two  ([34]).  The  numerical  experiments  for  the  computation  of 
harmonic  measure  in  two  dimensions  would  provide  experience  and  insights  for 
the  study  of  harmonic  measure  on  fractals  in  three  dimensions,  and  eventually 
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might  lead  to  the  understanding  of  harmonic  measure  on  porous  surfaces. 
Integral  Equations  On  Fractals 

In  growth  phenomena  such  as  crystallization,  electrodeposition,  viscous  finger¬ 
ing,  and  diffusion-limited  aggregation,  the  harmonic  measure  governs  the  growth 
of  the  fractal  surfaces,  in  addition  to  describing  the  distribution  of  growth  prob¬ 
abilities  ([56]). 

A  minor  modification  of  the  algorthm  of  this  thesis  can  be  used  to  study  the 
metric  properties  of  harmonic  measure  on  other  types  of  fractals. 

Integral  Equations  with  Other  Non-Oscillatory  Kernels 

It  is  easy  to  see  that  the  algorithm  of  this  thesis  does  not  substantially  depend 
on  the  fact  the  kernel  in  the  integral  equation  being  solved  satisfies  the  Laplace 
equation.  The  property  of  the  kernel  being  used  is  simply  its  smoothness  away 
from  the  diagonal,  and  the  fact  that  it  is  non-oscillatory.  Thus,  the  algorithm 
of  this  thesis  can  be  generalized  to  a  wide  class  of  integral  equations  both  in 
two  and  three  dimensions.  This  work  is  in  progress  and  will  be  reported  at  a 
later  date. 

Integral  Equations  On  Curves  in  R?  and  R? 

One  of  the  approaches  to  the  computation  of  electrostatic  fields  in  the  design  of 
chips  and  circuits  is  to  formulate  the  problems  as  integral  equations  on  curves 
either  in  FP  or  FP  with  a  free-space  Green’s  function  as  the  kernel  (see  [59]  and 
[47]).  The  algorithm  of  this  thesis  can  be  generalized  to  include  fast  algorithms 
for  problems  of  this  type.  The  new  algorithms  would  have  advantages  over  the 
conventional  methods  (such  as  moment  method,  and  finite  element  method) 
both  in  the  CPU  time  requirements,  and  the  accuracy  of  the  solution. 
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Appendix  A 
Measure  Theory 


In  this  appendix,  we  collect  relevant  facts  in  the  literature  about  Borel  integration, 
zind  harmonic  measure.  Facts  about  Borel  meaisure  can  be  found,  for  example,  in  [35], 
and  about  harmonic  measure  in  [21],  [15],  [31],  and  [55]. 


A.l  Borel  Measure  and  Integration 

Dehnition  A.l  (6-ring)  Let  TZ  be  a  non-empty  family  of  sets.  It  is  a  6-ring  if  the 
following  three  conditions  hold. 

1.  For  any  sets  A  ^IZ  and  B  ^  TZ, 


A[jB^n.  (A.l) 

2.  For  any  sets  A  ^IZ  and  B  ^7Z, 

A\BeTZ.  (A.2) 

3.  For  a  countable  sequence  of  sets  {A„}  C  TZ, 

PiAn^TZ.  (A.3) 
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In  other  words,  a  8-ring  is  a  nonempty  family  of  sets,  closed  under  difference,  finite 
union,  and  countable  intersection  operations. 

Definition  A. 2  (Measure)  Suppose  that  R  is  a  8-ring.  Then  a  set  function 

^  :  7?.  ^  [0,  oo) 

is  a  measure  if  for  any  disjoint  sequence  of  sets  {An}  C  R  with  \J  An  £R, 

(^-4) 

In  other  words,  a  measure  is  a  real,  non-negative,  and  countably  additive  function  on 
a  8 -ring. 

Definition  A. 3  (Borel  Ring  and  Borel  Measure)  Suppose  that  X  is  a  locally  compact 
Hausdorff  space,  and  R  is  the  8-ring  generated  by  the  family  of  compact  subsets  of  X. 
Then  the  8-ring  R  is  called  the  Borel  8-ring,  and  a  measure  on  R  is  called  a  Borel 
measure. 

Definition  A. 4  (Integrable  Function)  Suppose  that  X  is  a  locally  compact  Hausdorff 
space,  R  is  a  Borel  8-ring,  and  p  is  a  Borel  measure  on  R.  Then  a  function  f  :  X 
R^  is  p-integrable  if  there  exists  a  sequence  of  sets  {An}  C  R,  and  {a„}  €  R^  such 
that 

OO 

|an|/^(An)  <  OO,  (A.5) 

n=l 

and  for  any  x  Q.  X 

OO 

i/WI  <  Klx/in(2^),  (A.6) 

n=l 

where  xa„  is  the  characteristic  function  of  the  set  A„.  The  integration  of  function  f 
on  domain  X  with  respect  to  the  measure  p  is  defined  by 

OO 

fdp  =  lim  V  an,jp{An,j),  (A. 7) 

where  the  sum  satisfies  (A.5)  with  {onj}  €  R^  sets  {A„j}  C  R.  and  for  p  almost 
all  X, 


APPENDIX  A.  MEASURE  THEORY 


77 


and 


OO 


1  <  \f{x)l 

n=l 

(A.8) 

CO 

lim  Y,  (in,jXA„,,{x)  =  f{x). 

(A.9) 

A. 2  Harmonic  Measure 

Definition  A. 5  Suppose  that  fl  C  R^  is  a  region  with  F  as  its  boundary,  and  E  is  a 
closed  set  on  F.  Suppose  further  that  PziE)  is  the  solution  of  the  Dirichlet  problem 
for  D  with  the  boundary  value  1  on  E,  and  0  on  T\E ,  Then  the  solution  Pz{E)  is 
referred  to  as  the  harmonic  measure  of  E  at  the  point  z,  with  respect  to  Pi. 

Theorem  A.l  is  a  non-negative  measure,  and  Pz{U)  =  1- 

Theorem  A. 2  Suppose  that  p  and  q  are  two  points  in  fJ.  Then  there  exists  a  constant 
M  such  that 

pAE)  <  MpAE)  (A.io) 

for  any  measurable  set  E  CT. 


